人工蜂群 (Artificial Bee Colony, ABC)算法是土耳其学者 Karaboga D 为解决函数优化问题提于2005年提出的一种智能优化算法。其更快、收敛精度更高。人工蜂群算法因其性能良好, 受到人们广泛夹注, 已成为解决非线性连续优化问题的有效工具并已被成功应用到多个领域中。
对于蜜蜂这类群居昆虫, 单只蜜蜂的行为比较简单, 但由这种简单个体所组成的群体却表现出非常复杂且有条不紊的行为, 它们适应环境变化的能力较强, 并能较快聚集到该环境中的较好食物源处采蜜。通常在一个蜂巢中有三种类型的蜜蜂:蜂王 雄蜂和工蜂。 三种蜜蜂的职责各不相同:蜂王负责产卵; 雄蜂与蜂王交配笅殖后代:工蜂负责采蜜、筑巢、清洁、守卫等各项工作。其中, 以工蜂的行为最为复杂,一般地, 在一个蜂群中,大多数的工蜂都留在蜂巢内值 “内勤”,只有少数作为 “侦察员”四处寻找食物源, 这些 “侦察员”专门搜索新的食物源,一旦发现,它们立刻变成采集蜂, 飞回蜂巢跳上一支圆圈舞蹈或 8 字形舞蹈。蜜蜂这种传递信息的舞蹈被称为 “摇摆舞”,摇摆舞的持续时间暗 示食物源与蜂巢之间的距离, 而蜜蜂从一侧向另一侧摆动由其腹部构成的环形中轴相对于太阳的方向就是食物源的方向, 摇摆舞的剧烈程度反映食物源的质量,身上附着的花粉味道则反映食物源的种类。在这种信息的指引下, 整个蜂群趋于飞向质量最好的食物源采蜜。由此可见, 蜜蜂群体这种奇妙的受食方式不仅可以充分利用个体的功能特性, 而且也能最快地适应环境 (资源) 的变化,当已有食物源被耗尽,或者勘探到更优食物源(或环境变化产生)时, “侦察员” 可以通过传递信息的方式引导整个蜂群尽快飞向新的食物源。
人工蜂群算法正是模拟工蜂的这种受食行为提出的, 由三个基本部分组成, 包括食物源、 㢈佣蜂和未雇佣蜂, 并且定义两种行为,包括招募蜜蜂到食物源和放弄食物源。
为了更好地了解人工蜂群算法中招募蜜蜂到食物源和放弃食物源的行为, 以下图为例进行说明。
假设有两个已发现的食物源 $\mathrm{A}$ 和 $\mathrm{B}$ 。非雇佣蜂一开始没有任何周围食物源的信息, 侦察蜂鉴于内部和外部因素搜索周围的潜在食物源(过程 $\mathrm{S}$ ),若侦察蜂找到潜在食物源, 则利用自身能力记录下食物源相关信息, 返回蜂巢跳摇摆舞; 若招莫到蜜蜂采蜜, 则潜在食物源 变为新的食物源, 相应的侦察蜂转化为引领蜂(过程 $\mathrm{R}$ )。引领蜂从食物源(A 或 B)处携带花蜜返回蜂巢, 并将花蜜卸载到储存花蜜的位置, 在卸下花蜜后, 蜜蜂有以下三种选择:
(1) 放弃食物源成为末雇佣蜂 (过程 UF) ;
(2) 在返回到相同的食物源之前, 跳摇摆舞招募蜜蜂到该食物源处采蜜 (过程 EF1);
(3) 继续在该食物源采蜜而不招募任何蜜蜂 (过程 EF2)。 这里需要说明的是: 不是所有的蜜蜂都同时采蜜。
人工蜂群算法模拟实际蜜蜂采蜜机制处理函数优化问题, 该算法的基本思想是从某个随机产生的初始群体开始, 在适应度值较优的一半个体周围搜索, 采用一对一的竞争生存策略择优保留个体, 该操作称为引领蜂搜索。然后利用轮盘赌选择方式选择较优个体, 并在其周围进行贪婪搜索, 产生另一半较优个体, 这一过程称为跟随蜂搜索。将引领蜂和跟随蜂产生 的个体组成新的种群, 以避免种群多样性丧失, 并进行侦察蜂的类变异搜索,形成迭代种群。 该算法通过不断的迭代计算, 保留优质个体, 淘汰劣质个体, 向全局最优解靠近。
设置初始进化代数 $t=0$ ,在优化问题的可行解空间内随机产生满足约束条件 的 $\mathrm{NP}$ 个个体 $X$ 构成初始种群。 $$ X_i^0=\mathrm{Ub}+\operatorname{rand} \times(\mathrm{ub}-\mathrm{lb}), i=1,2, \cdots, \mathrm{NP} $$ 其中, $X_i^0$ 表示第 0 代种群中的第 $i$ 个个体; rand 为区间 $[0,1]$ 内的随机数。ub 为寻优参数的上边界, $1 \mathrm{~b}$ 为寻优参数的下边界。
种群中适应度值较小的一半个体构成引领蜂种群, 另一半个体构成跟随蜂种群。对于当前的第$t$代引领蜂种群中的一个目标个体$X_i^t$,随机选择个体$r_1\in [1,2,...,NP/2](r_1\neq i)$逐维进行交叉搜索,产生新个体$V$ $$ V(j)=X_i^t(j)+(2 \times \text { rand }-1) \times\left(X_i^t(j)-X_{r 1}^t(j)\right) $$ 其中, rand 为区间 $[0,1]$ 内的随机数。 如下图所示的是目标函数为二维时的交叉搜索示意图。
从图中可以看出, 目标引领蜂个体与随机选择引领蜂个体所形成的差分向量的方向和大小具有不确定性, 将此差分向量加至基向量上, 相当于在基向量上附加上一个规定范围内的随机扰动, 增加种群多样性。
与其他进化算法一样, 人工蜂群算法采用达尔文进化论 “优胜劣汰” 的思想来择优保留个体, 以保证算法不断向全局最优解进化。对新生成个体 $V$ 和目标个体 $X_i^t$ 进行适应度评价, 再将二者的适应度值进行比较,选择适应度值较优的个体进入引领蜂种群。 $$ X_i^{t+1}= \begin{cases}V & , f(V)<f\left(X_i^t\right) \\ X_i^t & , f\left(X_i^t\right) \geq f(V)\end{cases} $$ 其中 $f$ 为适应度函数。
跟随蜂搜索 以轮盘赌选择的方式在新的引领蜂种群中选择较优目标个佫 $X_k^{t+1}(k \in[1, \cdots \mathrm{NP} / 2])$, 与随机选择的个体进行搜索, 产生新个体 $X_k^{t+1}(k \in[\mathrm{NP} \mid$ $2+1, \cdots, \mathrm{NP})$, 并形成跟随蜂种群。 $$ P_i=\frac{\mathrm{fit}_i}{\sum_{i=1}^{\mathrm{NP} / 2} \mathrm{fit}_i} $$
人工蜂群算法中跟随蜂种群的搜索方式是其区别于其他进化算法的关键, 其本质上是择优选择个体进行贪婪搜索,是算法快速收敛的关键,但因为其搜索过程中引入随机信息, 所以不会过多减少种群多样性。
经过引领蜂种群和跟随蜂种群的搜索以后结合形成与初始种群规模大小相同的新种群, 为了避免随着种群进化种群多样性丧失过多, 人工蜂群算法模拟侦察蜂搜索潜在食物源的生理行为, 提出特有的侦察蜂搜索方式。假设某一个体连续 limit 代不变, 相应个体转换成侦察蜂, 按搜索产生新个体, 并与原个体按式 进行一对一比较, 保留适应度值较优的个体。
人工蜂群算法通过以上引领蜂种群、跟随蜂种群及侦察蜂种群的搜索, 使种群进化到下 一代并反复循环, 直到算法迭代次数 $t$ 达到预定的最大迭代次数 $G$ 时算法结束。
import numpy as np
import copy as copy
import numpy as np
from matplotlib import pyplot as plt
def initialization(pop,ub,lb,dim):
''' 种群初始化函数'''
'''
pop:为种群数量
dim:每个个体的维度
ub:每个维度的变量上边界,维度为[dim,1]
lb:为每个维度的变量下边界,维度为[dim,1]
X:为输出的种群,维度[pop,dim]
'''
X = np.zeros([pop,dim]) #声明空间
for i in range(pop):
for j in range(dim):
X[i,j]=(ub[j]-lb[j])*np.random.random()+lb[j] #生成[lb,ub]之间的随机数
return X
def BorderCheck(X,ub,lb,pop,dim):
'''边界检查函数'''
'''
dim:为每个个体数据的维度大小
X:为输入数据,维度为[pop,dim]
ub:为个体数据上边界,维度为[dim,1]
lb:为个体数据下边界,维度为[dim,1]
pop:为种群数量
'''
for i in range(pop):
for j in range(dim):
if X[i,j]>ub[j]:
X[i,j] = ub[j]
elif X[i,j]<lb[j]:
X[i,j] = lb[j]
return X
def CaculateFitness(X,fun):
'''计算种群的所有个体的适应度值'''
pop = X.shape[0]
fitness = np.zeros([pop, 1])
for i in range(pop):
fitness[i] = fun(X[i, :])
return fitness
def SortFitness(Fit):
'''适应度排序'''
'''
输入为适应度值
输出为排序后的适应度值,和索引
'''
fitness = np.sort(Fit, axis=0)
index = np.argsort(Fit, axis=0)
return fitness,index
def SortPosition(X,index):
'''根据适应度对位置进行排序'''
Xnew = np.zeros(X.shape)
for i in range(X.shape[0]):
Xnew[i,:] = X[index[i],:]
return Xnew
def RouletteWheelSelection(P):
'''轮盘赌策略'''
C = np.cumsum(P)#累加
r = np.random.random()*C[-1]#定义选择阈值,将随机概率与总和的乘积作为阈值
out = 0
#若大于或等于阈值,则输出当前索引,并将其作为结果,循环结束
for i in range(P.shape[0]):
if r<C[i]:
out = i
break
return out
def ABC(pop, dim, lb, ub, MaxIter, fun):
'''人工蜂群算法'''
'''
输入:
pop:为种群数量
dim:每个个体的维度
ub:为个体上边界信息,维度为[1,dim]
lb:为个体下边界信息,维度为[1,dim]
fun:为适应度函数接口
MaxIter:为最大迭代次数
输出:
GbestScore:最优解对应的适应度值
GbestPositon:最优解
Curve:迭代曲线
'''
L = round(0.6*dim*pop) #limit 参数
C = np.zeros([pop,1]) #计数器,用于与limit进行比较判定接下来的操作
nOnlooker=pop #引领蜂数量
X= initialization(pop,ub,lb,dim) # 初始化种群
fitness = CaculateFitness(X, fun) # 计算适应度值
fitness, sortIndex = SortFitness(fitness) # 对适应度值排序
X = SortPosition(X, sortIndex) # 种群排序
GbestScore = copy.copy(fitness[0])#记录最优适应度值
GbestPositon = np.zeros([1,dim])
GbestPositon[0,:] = copy.copy(X[0, :])#记录最优位置
Curve = np.zeros([MaxIter, 1])
Xnew = np.zeros([pop,dim])
fitnessNew = copy.copy(fitness)
for t in range(MaxIter):
'''引领蜂搜索'''
for i in range(pop):
k = np.random.randint(pop)#随机选择一个个体
while(k==i):#当k=i时,再次随机选择,直到k不等于i
k = np.random.randint(pop)
phi = (2*np.random.random([1,dim]) - 1)
Xnew[i,:] = X[i,:]+phi*(X[i,:]-X[k,:]) #公式(2.2)位置更新
Xnew=BorderCheck(Xnew,ub,lb,pop,dim) #边界检查
fitnessNew = CaculateFitness(Xnew, fun) # 计算适应度值
for i in range(pop):
if fitnessNew[i]<fitness[i]:#如果适应度值更优,替换原始位置
X[i,:]= copy.copy(Xnew[i,:])
fitness[i] = copy.copy(fitnessNew[i])
else:
C[i] = C[i]+1 #如果位置没有更新,累加器+1
#计算选择适应度权重
F = np.zeros([pop,1])
MeanCost = np.mean(fitness)
for i in range(pop):
F[i]=np.exp(-fitness[i]/MeanCost)
P=F/sum(F) #式(2.4)
'''侦察蜂搜索'''
for m in range(nOnlooker):
i=RouletteWheelSelection(P)#轮盘赌测量选择个体
k = np.random.randint(pop)#随机选择个体
while(k==i):
k = np.random.randint(pop)
phi = (2*np.random.random([1,dim]) - 1)
Xnew[i,:] = X[i,:]+phi*(X[i,:]-X[k,:])#位置更新
Xnew=BorderCheck(Xnew,ub,lb,pop,dim)#边界检查
fitnessNew = CaculateFitness(Xnew, fun) # 计算适应度值
for i in range(pop):
if fitnessNew[i]<fitness[i]:#如果适应度值更优,替换原始位置
X[i,:]= copy.copy(Xnew[i,:])
fitness[i] = copy.copy(fitnessNew[i])
else:
C[i] = C[i]+1 #如果位置没有更新,累加器+1
'''判断limit条件,并进行更新'''
for i in range(pop):
if C[i]>=L:
for j in range(dim):
X[i, j] = np.random.random() * (ub[j] - lb[j]) + lb[j]
C[i] = 0
fitness = CaculateFitness(X, fun) # 计算适应度值
fitness, sortIndex = SortFitness(fitness) # 对适应度值排序
X = SortPosition(X, sortIndex) # 种群排序
if fitness[0] <= GbestScore: # 更新全局最优
GbestScore = copy.copy(fitness[0])
GbestPositon[0,:] = copy.copy(X[0, :])
Curve[t] = GbestScore
return GbestScore, GbestPositon, Curve
'''人工蜂群优化算法求解压力容器设计'''
'''适应度函数'''
def fun(X):
x1 = X[0] #Ts
x2 = X[1] #Th
x3 = X[2] #R
x4 = X[3] #L
#约束条件判断
g1 = -x1+0.0193*x3
g2 = -x2+0.00954*x3
g3 = -np.math.pi*x3**2-4*np.math.pi*x3**3/3+1296000
g4 = x4-240
if g1<=0 and g2<=0 and g3<=0 and g4<=0:
#如果满足约束条件则计算适应度值
fitness = 0.6224*x1*x3*x4 + 1.7781*x2*x3**2 + 3.1661*x1**2*x4 + 19.84*x1**2*x3
else:
#如果不满足约束条件,则设置适应度值为很大的一个惩罚数
fitness = 10E32
return fitness
'''主函数 '''
#设置参数
pop = 50 #种群数量
MaxIter = 500 #最大迭代次数
dim = 4 #维度
lb = np.array([0,0,10,10]) #下边界
ub = np.array([100,100,100,100])#上边界
#适应度函数选择
fobj = fun
GbestScore,GbestPositon,Curve = ABC(pop,dim,lb,ub,MaxIter,fobj)
print('最优适应度值:',GbestScore)
print('最优解[Ts,Th,R,L]:',GbestPositon)
#绘制适应度曲线
plt.figure(1)
plt.plot(Curve,'r-',linewidth=2)
plt.xlabel('Iteration',fontsize='medium')
plt.ylabel("Fitness",fontsize='medium')
plt.grid()
plt.title('ABC',fontsize='large')
plt.show()
参考资料:
import numpy as np
import copy as copy
import numpy as np
from matplotlib import pyplot as plt
def initialization(pop,ub,lb,dim):
''' 种群初始化函数'''
'''
pop:为种群数量
dim:每个个体的维度
ub:每个维度的变量上边界,维度为[dim,1]
lb:为每个维度的变量下边界,维度为[dim,1]
X:为输出的种群,维度[pop,dim]
'''
X = np.zeros([pop,dim]) #声明空间
for i in range(pop):
for j in range(dim):
X[i,j]=(ub[j]-lb[j])*np.random.random()+lb[j] #生成[lb,ub]之间的随机数
return X
def BorderCheck(X,ub,lb,pop,dim):
'''边界检查函数'''
'''
dim:为每个个体数据的维度大小
X:为输入数据,维度为[pop,dim]
ub:为个体数据上边界,维度为[dim,1]
lb:为个体数据下边界,维度为[dim,1]
pop:为种群数量
'''
for i in range(pop):
for j in range(dim):
if X[i,j]>ub[j]:
X[i,j] = ub[j]
elif X[i,j]<lb[j]:
X[i,j] = lb[j]
return X
def CaculateFitness(X,fun):
'''计算种群的所有个体的适应度值'''
pop = X.shape[0]
fitness = np.zeros([pop, 1])
for i in range(pop):
fitness[i] = fun(X[i, :])
return fitness
def SortFitness(Fit):
'''适应度排序'''
'''
输入为适应度值
输出为排序后的适应度值,和索引
'''
fitness = np.sort(Fit, axis=0)
index = np.argsort(Fit, axis=0)
return fitness,index
def SortPosition(X,index):
'''根据适应度对位置进行排序'''
Xnew = np.zeros(X.shape)
for i in range(X.shape[0]):
Xnew[i,:] = X[index[i],:]
return Xnew
def RouletteWheelSelection(P):
'''轮盘赌策略'''
C = np.cumsum(P)#累加
r = np.random.random()*C[-1]#定义选择阈值,将随机概率与总和的乘积作为阈值
out = 0
#若大于或等于阈值,则输出当前索引,并将其作为结果,循环结束
for i in range(P.shape[0]):
if r<C[i]:
out = i
break
return out
def ABC(pop, dim, lb, ub, MaxIter, fun):
'''人工蜂群算法'''
'''
输入:
pop:为种群数量
dim:每个个体的维度
ub:为个体上边界信息,维度为[1,dim]
lb:为个体下边界信息,维度为[1,dim]
fun:为适应度函数接口
MaxIter:为最大迭代次数
输出:
GbestScore:最优解对应的适应度值
GbestPositon:最优解
Curve:迭代曲线
'''
L = round(0.6*dim*pop) #limit 参数
C = np.zeros([pop,1]) #计数器,用于与limit进行比较判定接下来的操作
nOnlooker=pop #引领蜂数量
X= initialization(pop,ub,lb,dim) # 初始化种群
fitness = CaculateFitness(X, fun) # 计算适应度值
fitness, sortIndex = SortFitness(fitness) # 对适应度值排序
X = SortPosition(X, sortIndex) # 种群排序
GbestScore = copy.copy(fitness[0])#记录最优适应度值
GbestPositon = np.zeros([1,dim])
GbestPositon[0,:] = copy.copy(X[0, :])#记录最优位置
Curve = np.zeros([MaxIter, 1])
Xnew = np.zeros([pop,dim])
fitnessNew = copy.copy(fitness)
for t in range(MaxIter):
'''引领蜂搜索'''
for i in range(pop):
k = np.random.randint(pop)#随机选择一个个体
while(k==i):#当k=i时,再次随机选择,直到k不等于i
k = np.random.randint(pop)
phi = (2*np.random.random([1,dim]) - 1)
Xnew[i,:] = X[i,:]+phi*(X[i,:]-X[k,:]) #公式(2.2)位置更新
Xnew=BorderCheck(Xnew,ub,lb,pop,dim) #边界检查
fitnessNew = CaculateFitness(Xnew, fun) # 计算适应度值
for i in range(pop):
if fitnessNew[i]<fitness[i]:#如果适应度值更优,替换原始位置
X[i,:]= copy.copy(Xnew[i,:])
fitness[i] = copy.copy(fitnessNew[i])
else:
C[i] = C[i]+1 #如果位置没有更新,累加器+1
#计算选择适应度权重
F = np.zeros([pop,1])
MeanCost = np.mean(fitness)
for i in range(pop):
F[i]=np.exp(-fitness[i]/MeanCost)
P=F/sum(F) #式(2.4)
'''侦察蜂搜索'''
for m in range(nOnlooker):
i=RouletteWheelSelection(P)#轮盘赌测量选择个体
k = np.random.randint(pop)#随机选择个体
while(k==i):
k = np.random.randint(pop)
phi = (2*np.random.random([1,dim]) - 1)
Xnew[i,:] = X[i,:]+phi*(X[i,:]-X[k,:])#位置更新
Xnew=BorderCheck(Xnew,ub,lb,pop,dim)#边界检查
fitnessNew = CaculateFitness(Xnew, fun) # 计算适应度值
for i in range(pop):
if fitnessNew[i]<fitness[i]:#如果适应度值更优,替换原始位置
X[i,:]= copy.copy(Xnew[i,:])
fitness[i] = copy.copy(fitnessNew[i])
else:
C[i] = C[i]+1 #如果位置没有更新,累加器+1
'''判断limit条件,并进行更新'''
for i in range(pop):
if C[i]>=L:
for j in range(dim):
X[i, j] = np.random.random() * (ub[j] - lb[j]) + lb[j]
C[i] = 0
fitness = CaculateFitness(X, fun) # 计算适应度值
fitness, sortIndex = SortFitness(fitness) # 对适应度值排序
X = SortPosition(X, sortIndex) # 种群排序
if fitness[0] <= GbestScore: # 更新全局最优
GbestScore = copy.copy(fitness[0])
GbestPositon[0,:] = copy.copy(X[0, :])
Curve[t] = GbestScore
return GbestScore, GbestPositon, Curve
'''人工蜂群优化算法求解压力容器设计'''
'''适应度函数'''
def fun(X):
x1 = X[0] #Ts
x2 = X[1] #Th
x3 = X[2] #R
x4 = X[3] #L
#约束条件判断
g1 = -x1+0.0193*x3
g2 = -x2+0.00954*x3
g3 = -np.math.pi*x3**2-4*np.math.pi*x3**3/3+1296000
g4 = x4-240
if g1<=0 and g2<=0 and g3<=0 and g4<=0:
#如果满足约束条件则计算适应度值
fitness = 0.6224*x1*x3*x4 + 1.7781*x2*x3**2 + 3.1661*x1**2*x4 + 19.84*x1**2*x3
else:
#如果不满足约束条件,则设置适应度值为很大的一个惩罚数
fitness = 10E32
return fitness
'''主函数 '''
#设置参数
pop = 50 #种群数量
MaxIter = 500 #最大迭代次数
dim = 4 #维度
lb = np.array([0,0,10,10]) #下边界
ub = np.array([100,100,100,100])#上边界
#适应度函数选择
fobj = fun
GbestScore,GbestPositon,Curve = ABC(pop,dim,lb,ub,MaxIter,fobj)
print('最优适应度值:',GbestScore)
print('最优解[Ts,Th,R,L]:',GbestPositon)
#绘制适应度曲线
plt.figure(1)
plt.plot(Curve,'r-',linewidth=2)
plt.xlabel('Iteration',fontsize='medium')
plt.ylabel("Fitness",fontsize='medium')
plt.grid()
plt.title('ABC',fontsize='large')
plt.show()
最优适应度值: [8050.93401837] 最优解[Ts,Th,R,L]: [[ 1.30055305 0.64286315 67.38602337 10.0000763 ]]